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Abstract 

We  present  a new  approach  to  the  construction  of  biorthogonal  wavelet  transforms  us- 
ing polynomial  splines.  The  construction  is  performed  in  a “lifting”  manner  and  we  use 
interpolatory,  as  well  as  local  quasi-interpolatory  and  smoothing  splines  as  predicting 
aggregates  in  this  scheme.  The  transforms  contain  some  scalar  control  parameters  which 
enable  their  flexible  tuning  in  either  time  or  frequency  domains.  The  transforms  are 
implemented  in  a fast  way.  They  demonstrated  efficiency  in  application  to  image  com- 
pression. 


1 Introduction 

Until  recently,  two  methods  have  been  used  for  the  construction  of  wavelet  schemes 
using  splines.  One  is  to  construct  orthogonal  and  semi-orthogonal  wavelets  in  the  spline 
spaces  (Battle-Lemarie  [2,  7],  Chui-Wang  [6],  Unser-Aldroubi-Eden  [12]).  Another  way 
was  introduced  by  Cohen,  Daubechies  and  Feauveau  [3]  who  constructed  symmetric 
compactly  supported  spline  wavelets  whose  duals,  remaining  compactly  supported  and 
symmetric,  do  not  belong  to  a spline  space.  However,  since  the  introduction  of  the  lifting 
scheme  for  the  design  of  wavelet  transforms  [11],  a new  way  was  opened  to  use  splines  as 
a tool  for  devising  a full  discrete  scheme  of  wavelet  transforms.  Namely,  various  splines 
can  be  employed  as  predicting  aggregates  in  lifting  constructions. 

2 Lifting  scheme  of  biorthogonal  wavelet  transform 

The  sequences  which  belong  to  the  space  Zj,  we  call  the  discrete-time  sig- 
nals. The  ^-transform  of  a signal  {a(A;)}  is  defined  as  follows:  a(z)  = o.{k). 

Throughout  the  paper  we  assume  that  z = e’".  We  introduce  a family  of  biorthogonal 
wavelet-type  transforms  that  operate  on  the  signal  x = which  we  construct 

through  lifting  steps. 

The  lifting  scheme  for  the  wavelet  transform  of  a signal  can  be  implemented  in  primal 
or  dual  modes.  For  brevity  we  consider  only  the  primal  mode. 

Decomposition  Generally,  the  primal  lifting  scheme  for  decomposition  of  signals  con- 
sists of  three  steps:  1.  Split.  2.  Predict.  3.  Update  or  lifting. 

Split  - We  split  the  array  x into  even  and  odd  sub-arrays: 

= {ei(fc)  = x{2k)},  di  = {di{k)  = x{2k  -f  1)},  fc  G Z. 
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Predict  - We  use  the  even  array  ei  to  predict  the  odd  array  di  and  redefine  the 
array  di  as  the  difference  between  the  existing  array  and  the  predicted  one.  To  be  specific, 
we  apply  some  filter  with  transfer  function  zU{z)  to  the  sequence  ej  and  predict  the 
function  di(z^)  which  is  the  2:^— transform  of  di.  The  2^— transform  of  the  new  d— array 
is  defined  as  follows: 

di{z^)  = di{z^)  — zU{z)ei{z^).  (2.1) 

Prom  now  on  the  superscript  u means  an  update  operation  of  the  array.  Obviously,  the 
prediction  zU{z)ei{z'^)  should  approximate  di(z^)  well. 

Lifting  - We  update  the  even  array  using  the  new  odd  array; 

eUz^)^ei{z^)  + Piz)^-^d^{z^)-  (2-2) 

Generally,  the  goal  of  this  step  is  to  eliminate  aliasing  which  appears  while  downsampling 
the  original  signal  x into  ei . Further  on  we  will  discuss  how  to  achieve  this  effect  by  a 
proper  choice  of  the  filter  /3. 

Reconstruction  The  reconstruction  of  the  signal  x from  the  arrays  e"  and  d“  is 
implemented  in  reverse  order:  1.  Undo  Lifting.  2.  Undo  Predict.  3.  Unsplit. 

Undo  Lifting  - We  restore  the  even  array;  ei(z^)  = e'y(z^)  - /3{z)z~^  di{z'^). 
Undo  Predict  - We  restore  the  odd  array:  d\{z^)  = + zU{z)ei{z^). 

Unsplit  - The  last  step  represents  the  standard  restoration  of  the  signal  from  its 
even  and  odd  components.  In  the  z-domain  this  is  x{z)  = ei(z^)  + z~^di{z'^). 

The  lifting  scheme  presented  above,  yields  an  efficient  algorithm  for  the  implementa- 
tion of  the  forward  and  backward  transform  of  x > e“  U d“.  These  operations  can  be 
interpreted  as  a transformation  of  the  signal  by  a filter  bank  that  possesses  the  perfect 
reconstruction  properties  and  it  is  associated  with  the  biorthogonal  pairs  of  bases  in  the 
space  of  discrete-time  signals.  These  basis  signals  are  synthesis  and  analysis  wavelets. 
Further  steps  of  the  transform  are  implemented  in  an  iterative  way  by  the  same  lifting 
operations. 

3 Polynomial  splines 

We  will  construct  polynomial  splines  of  various  kinds  using  the  even  subarray  of  a signal, 
calculate  their  values  in  the  midpoints  between  nodes  and  use  these  values  for  prediction 
of  the  odd  array.  In  this  section  we  discuss  some  properties  of  such  splines  and  derive 
the  corresponding  filters  U. 

3.1  B— splines 

The  central  B— spline  of  first  order  on  the  grid  {kh}  is  defined  as  follows: 

= / ifx€  [-/i/2,fi/2], 

^ ^ (0  elsewhere. 

The  central  B— spline  of  order  p is  the  convolution  (x)  = (x)  * M\  (x)  p > 2. 

Note  that  the  B— spline  of  orderp  is  supported  at  the  interval  {—ph/2,ph/2).  It  is  positive 
within  its  support  and  symmetric  around  zero.  The  nodes  of  B— splines  of  even  orders 
are  located  at  points  {kfi}  and  of  odd  orders  at  points  {h{k  + 1/2)},  fc  € Z.  It  is  readily 
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verified  that  hMf^{hx)  = Mp{x),  where  M'P{x)  Mf(a;).  Let 

uP  :=  {hMlihk)  = M^{k)},  and  wP  :=  {hMl{h{k  + 1/2))  = M^{k  + 1/2)}  , fc  e Z. 

(3.1) 

Due  to  the  compact  support  of  5— splines,  these  sequences  are  finite.  We  will  use  for  our 
constructions  only  splines  of  odd  orders  p = 2r  - 1.  In  Table  1 we  present  the  sequences 
for  initial  values  r which  are  of  practical  importance. 


k 

-3 

-2 

-1 

0 

1 

2 

3 

X 8 

0 

0 

1 

6 

1 

0 

0 

u®  X 384 

0 

1 

76 

230 

76 

1 

0 

X 2 

0 

0 

1 

1 

0 

0 

0 

w®  X 24 

0 

1 

11 

11 

1 

0 

0 

Tab.  1.  Values  of  the  sequences  and  w^. 


We  need  the  transforms  of  the  sequences  uP  and  wP  : 

OO  OC 

uP(z2);=  Y,  z-^^u^{k),  wP{z^):=  Y 

/c=  — OO  k=—oo 

These  functions  are  Laurent  polynomials,  and  are  called  the  Euler-Frobenius  polynomials 

[10]. 

Proposition  3.1.  ([9])  On  the  circle  z = e‘“  the  Laurent  polynomials  u'’’{z^)  are 
strictly  positive.  Their  roots  are  all  simple  and  negative.  Each  root  C can  be  paired  with 
a dual  root  6 such  that  (^6  = 1.  Thus,  ifp  = 2r  + 1 is  odd,  then  uP(z^)  can  be  represented 
as  follows: 

^ 1 

==  n + TnZ^){l  + 0 < 7„  < 1.  (3.2) 

n=l  O'" 


We  denote 


Uf{z) 


iwPjz^) 

uP{z'^) 


Proposition  3.2  The  rational  functions  Uf{z)  are  real-valued  and  Uf{—z) 
Ifp  = 2r  + 1 is  odd  then 


l-Uf{z) 


{a-2r+%{a)  i-a-2y+%{-a) 

l + = up-(^)— ■ 


(3.3) 
= -Uf{z). 

(3.4) 


where  a :=  z -h  z ^ and  ^r(o)  is  a polynomial  of  degree  r — 1. 


3.2  Inter polatory  splines 

The  shifts  of  15— splines  form  a basis  in  the  space  of  splines  of  order  p on  the  grid  kh. 
Namely,  any  spline  Sf^  £ S]]  has  the  following  representation: 

Sl{x)  =hY^{l)K{^-lh).  (3.5) 

l 


Splines  and  wavelets 


317 


Let  q :=  {9(0};  q{z^)  be  the  transform  of  q.  We  introduce  also  the  se- 

quences :=  h{S^{hk)  = -Sf  (A;)}  and  mP  ;=  {S^{h{k  + 1/2))  — 5f(fc  -f  1/2)}  of  values 
of  the  spline  on  the  grid  points  and  on  the  midpoints.  Let  and  mP{z^)  be  the 

corresponding  a:^-transforms.  We  have 

S\{k)  = ^ q{l)  Ml{k  - 1),  and  5?  ^ q(l)  ^ -f  . (3.6) 

Respectively,  sP(^^)  = g(2:^)u(.z^),  and  wP{z^)  = q{z^)w{z^). 

Prom  these  formulae  we  can  derive  expression  for  the  coefficients  of  a spline  which 
interpolates  a given  sequence  e :=  {e(fc)}  at  grid  points; 

hSl(hk)  = e{k),  keZ,^  q{z^)vP{z)  = e(z^)  q(z^)  = • (3-7) 

The  transform  of  the  sequence  is: 

mP{z^)  = q{z‘^)w‘‘ {z"^)  = zUf{z)e{z^).  (3.8) 

Our  further  construction  exploits  the  super-convergence  property  of  the  interpolatory 
splines  of  odd  orders  (even  degrees). 

Theorem  3.3.  ([13])  Let  a function  f € L^(— oo,oo)  havep+1  continuous  derivatives 
and  let  interpolate  f on  the  grid  {kh}.  Denote  fk  = f{.{k  + l/2)h).  Then  in  the 

case  of  odd  p = 2r  + 1,  the  following  asymptotic  relation  holds. 

Sl{h{k+l/2))  = f,-h^r+2f(2r+2)^f^^j._^y2)){2r+l^^^^^^  +o{h^^+^ 

(3.9) 

where  bs{x)  is  the  Bernoulli  polynomial  of  degree  s. 

Recall,  that  in  general  the  interpolatory  spline  of  order  2r  -f  1 approximates  the 
function  / with  accuracy  of  Therefore,  we  may  claim  that  {(A;  -f  l/2)/i}  are 

points  of  super-convergence  of  the  spline  5^.  Note,  that  the  spline  of  order  2r  -|- 1,  which 
interpolates  the  values  of  a polynomial  of  degree  2r,  coincides  with  this  polynomial. 
However,  the  spline  of  order  2r  + 1 which  interpolates  the  values  of  a polynomial  of 
degree  2r  + 1 on  the  grid  {kh}  restores  the  values  of  this  polynomial  at  the  mid-points 
{(A;-M/2)/i}.  This  property  will  result  in  the  vanishing  moments  property  of  the  wavelets 
to  be  constructed  later. 

3.3  Quasi-inter polatory  splines 

We  can  see  from  (3.7)  and  (3.8)  that  in  order  to  find  values  at  the  midpoints  of  the  spline 
interpolating  the  signal  e,  the  signal  has  to  be  filtered  with  the  filter  whose  transfer 
function  is  zUf{z).  This  filter  has  infinite  impulse  response  (HR).  However,  the  property 
of  super-convergence  at  the  midpoints  is  not  an  exclusive  attribute  of  the  interpolatory 
splines.  It  is  also  inherent  to  the  so  called  local  quasi-interpolatory  splines  of  odd  orders, 
which  can  be  constructed  using  finite  impulse  response  (FIR)  filtering. 

Definition  3.4  Let  the  function  f have  p continuous  derivatives  and  f :=  {fk  = 
f{hk)},  A:  € Z.  The  spline  of  order  p given  by  (3.5)  is  said  to  be  the  local 
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quasi-interpolatory  spline  if  the  array  q of  its  coefficients  is  derived  by  FIR  filtering  the 
array  of  samples  f 

(3.10) 

where  r{z'^)  is  a Laurent  polynomial,  and  the  difference  \f{x)  — 5'^(a;)|  = 0{f^^^hP).  If 
f is  a polynomial  of  degree  p — I,  then  the  spline  Sf,{x)  = f(x). 

If  wP  is  the  sequence  defined  in  (3.1)  then  the  midpoint  values  m’’  are  produced  by 
the  following  FIR  filtering  of  the  array  of  samples  f:  mP(2;^)  = zUP{z)i{z^),  U^[z)  := 
z~^T{z'^)w^{z^).  Explicit  formulas  for  the  construction  of  quasi-interpolatory  splines  as 
well  as  the  estimations  of  the  differences  were  established  in  [13].  In  the  present  work 
we  are  interested  in  splines  of  odd  orders  p = 2r  -f  1.  There  are  many  FIR  filters  which 
generate  quasi-interpolatory  splines  but  only  one  filter  of  minimal  length  2r  -|- 1 for  each 
order  p = 2r  + 1.  Let  A(^)  :=  z~^  — 2-1-2:^. 

Theorem  3.5  A quasi-interpolatory  spline  of  order  p = 2r  -|-  1 can  be  produced  by 
filtering  (3.10)  with  filters  F of  length  no  less  than  2r  1.  There  exists  a unique  filter 
of  length  2r-|- 1 which  produces  the  minimal  quasi-interpolatory  spline  . Its 

transfer  function  is: 

r /„  . 2r-|-l  oo 

Tl,{z'^)  = l + ^Pl\\z),  =Y^{-l)^0lt^K  (3.11) 

fc=l  ' fc=0 

If  the  function  f has  2r  + 3 derivatives  then  the  following  asymptotic  relations  hold 
for  the  midpoint  values  of  the  minimal  quasi-interpolatory  spline  of  odd  order: 

Sl'-+\h{k  + 1/2))  = f{h{k  + 1/2))  -P  fi2'-+2/(2r+2)(^(^  ^ 1/2))A’-  + 0(/'2^+®Vl2’-+!’), 

(2r-f  l)&2r+2(0)  p. 

where  bg(x)  is  the  Bernoulli  polynomial  of  degree  s. 

This  implies  that  the  super-convergence  property  is  similar  to  that  of  the  interpol- 
atory  splines.  The  asymptotic  representation  (3.12)  provides  tools  for  custom  design  of 
predicting  splines  retaining  or  even  enhancing  the  approximation  accuracy  of  the  min- 
imal spline  at  the  midpoints. 

Proposition  3.6  If  the  coefficients  of  the  spline  order  2r  -j-1  are 

derived  as  in  (3.10)  using  the  filter  F^  of  length  2r  -t-  3,  with  the  transfer  function 
Fp(z^)  = V(„{z'^)  -f-  pA^'*'^(2:),  then  the  spline  restores  polynomials  of  degree  2r  -I-  1 at 
the  midpoints  between  nodes,  for  any  real  value  p.  However,  if  p = —A^  then  the  spline 
restores  polynomials  of  degree  2r  4-  3. 

If  the  parameter  p is  chosen  such  that  p = (— l)^jp|  then  the  spline  possesses 

the  smoothing  property  [14]. 
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3.4  Examples 
3.4.1  Quadratic  splines 

Interpolatory  spline  Let  a = z~^  + z.  Then 

4a 


U}{z) 


z'^  + 6 + z~^  ’ 
Minimal  spline  The  filters  are 

. 1 


rU^‘) 

and  l-[/^(z) 

Extended  spline 


8 


Mz),  Ui(z)  ^ 


-z-^  + 9z~^  +9z 


-1 


16 


{a-2f{z~^+A  + z) 


16 


T\{z)  = Vl{z-^)  + -\\z),Ul{z) 


3z-^  - 25z-3  + 150Z-1  + 150z  - 25z^  + 3x® 


256 


and  1 - t/g  (z)  = 


(a  - 2)3 (3^-2  + 18z-i  + 38  + 18z  + 3z2) 


256 


Remark  3.7  In  [5]  Donoho  presented  a scheme  where  an  odd  sample  is  predicted  by 
the  value  in  the  central  point  of  the  polynomial  of  odd  degree  which  interpolates  adjacent 
even  samples.  One  can  observe  that  our  filter  coincides  with  the  filter  derived  by 
Donoho ’s  scheme  using  the  cubic  interpolatory  polynomial.  The  filter  Ul  coincides  with 
the  filter  derived  using  the  interpolatory  polynomial  of  fifth  degree.  On  the  other  hand, 
the  filter  Ul  is  closely  related  to  the  commonly  used  Butterworth  filter  [8].  Namely,  in 
this  case  the  filter  transfer  functions  {l  + Ul(z))/2,  :=  {l-Ul{z))/2 

coincide  with  magnitude  squared  of  the  transfer  functions  of  the  discrete-time  low-pass 
and  high-pass  half-band  Butterworth  filters  of  order  4>  respectively. 

3.4.2  Splines  of  fifth  order  (fourth  degree) 

Interpolatory  spline 

(a -2)3  (a -10) 


2 , 16(^3  + llz  + lU-i  + z-3) 


U({z) 


, 1 - Uf{z) 


+ 76^2  + 230  + 76^-2  + ' 


+ 76^2  + 230  + 76z-2  + ^-4 
Minimal  spline  The  filter  is 

47(z-’’  + z"^)  + 89{z-3  + z^)  - 2277(z-3  + z^)  + 15965a 


Ui{z) 


27648 


4 Wavelet  transforms  using  spline  filters 

4.1  Choosing  the  filters  for  the  lifting  step 

In  the  previous  section  we  presented  a family  of  filters  U for  the  predicting  step  which 
were  originated  from  splines  of  various  types.  But,  as  it  is  seen  from  (2.2),  to  accomplish 
the  transform  we  have  to  define  the  filter  j3.  There  is  a remarkable  freedom  in  the  choice  of 
these  filters.  The  only  requirement  needed  to  guarantee  a perfect  reconstruction  property 
of  the  transform  is  that  P{—z)  = ^{z).  In  order  to  make  synthesis  and  analysis  filters 
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similar  in  their  properties,  we  choose  fi{z)  — U[z)/2,  where  U means  one  of  filters  U 
presented  above.  In  particular,  U may  coincide  with  the  filter  U which  was  used  for  the 
prediction. 

We  say  that  a wavelet  'ip  has  m vanishing  moments  if  the  following  relations  hold: 

Proposition  4.1  Suppose  the  filtersU{z)  and  0{z)  = tj{z)/2  are  used  for  the  predicting 
and  lifting  steps,  respectively.  If  1 — U{z)  contains  the  factor  [z  -2-\-\jzY  then  the 
high-frequency  analysis  wavelets  t/A  have  2r  vanishing  moments.  If,  in  addition  1 — U{z) 
contains  the  factor  {z-2-\-l/zY  then  the  synthesis  wavelet  has  2q  vanishing  moments, 
where  q = mm{p,r}. 

4.2  Implementation  of  the  transforms 

Suppose,  we  have  chosen  the  filter  P = U/2.  The  functions  zU(z)  and  zU{z)  depend  on 
z'^  and  we  write  F{z'^)  zU(z)  and  F(z'^)  zU(z).  Then  the  decomposition  procedure 
is  (see  (2.1),  (2.2)): 

d^^{z)  = di{z)  - F{z)ei{z),  e^{z)  = ei{z) ^F{z)  d'p(z).  (4.1) 

Equation  (4.1)  means  that  in  order  to  obtain  the  detail  array  d)' , we  must  process  the 
even  array  ei  with  the  filter  F,  with  transfer  function  F{z),  and  extract  the  filtered  array 
from  the  odd  array  di.  In  order  to  obtain  the  smoothed  array  e",  we  must  process  the 
detail  array  d“  with  the  filter  4>  that  has  the  transfer  function  4>(2)  = z~^F{z)/2  and 
add  the  filtered  array  to  the  even  array  ei.  But  the  filter  4>  differs  from  Fr/2  only  by 
one-sample  delay  and  it  operates  similarly.  Thus,  both  operations  of  the  decomposition 
are,  in  principle,  identical.  For  the  reconstruction  the  same  operation  is  conducted  in 
reverse  order. 

Therefore,  it  is  sufficient  to  outline  the  implementation  of  the  filtering  with  the  func- 
tion F{z). 

Implementation  of  FIR  filters  originating  from  local  splines  is  straightforward  and, 
therefore  we  only  make  a few  remarks  on  HR  filters  originating  from  interpolatory 
splines.  A detailed  description  can  be  found  in  [1].  Equations  (3.2)  and  (3.3)  imply 
that,  while  the  interpolatory  spline  of  order  2r  + 1 is  used,  the  transfer  function  F{z)  = 
P{z)/YYn=\  :^(1  + 7n-2J)(l  + 7n-2^~^)>  where  P{z)  is  the  Laurent  polynomial.  It  means 
that  the  HR  filter  F can  be  split  into  a cascade  consisting  of  a FIR  filter  with  the 
transfer  function  P(^),  r elementary  causal  recursive  filters  denoted  by  R{n],  and  r ele- 
mentary anti-causal  recursive  filters,  denoted  by  )?(n).  The  causal  and  anti-causal  filters 
operate  as  follows: 

y = P(njx  <t=»  y{l)  = x{l)  -f  7„y(l  - 1),  y = k{n)x.  y{l)  = x{l)  + 7„y(f  -I- 1). 

Example  4.2  (Example  of  recursive  filter)  We  present  IIR  filters  derived  from  the 
interpolatory  splines  of  third  order. 
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Let  = 3 — 2\/2  i 


0.172.  Then 


1 + z 


The  filter  can  be  implemented  with  the  following  cascade; 


xo(k)  = 47i(a;(A:)+a:(fc  + l)),  xi(k)  = xo(k)-'Y}xi(k- 1),  y{k)  = xi{k)-'fly{k  + l). 
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